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Abstract 

A library of C functions yielding exact solutions of potential and flux influences due to 
uniform surface distribution of singularities on flat triangular and rectangular elements has 
been developed. This library, ISLES, has been used to develop the neBEM solver that is both 
precise and fast in solving a wide range of problems of scientific and technological interest. 
Here we present the exact expressions proposed for computing the influence of uniform 
singularity distributions on triangular elements and illustrate their accuracy. We also present 
a study concerning the time taken to evaluate these long and complicated expressions vis 
a vis that spent in carrying out simple quadratures. Finally, we solve a classic benchmark 
problem in electrostatics, namely, estimation of the capacitance of a unit square plate raised 
to unit volt. For this problem, we present the estimated values of capacitance and compare 
them successfully with some of the most accurate results available in the literature. In 
addition, we present the variation of the charge density close to the corner of the plate for 
various degrees of discretization. The variations are found to be smooth and converging. 
This is in clear contrast to the criticism commonly leveled against usual BEM solvers. 

Keywords: Boundary element method, triangular element, potential, flux, unit square 
plate, charge density, capacitance. 

1 Introduction 

One of the elegant methods for solving the Laplace / Poisson equations (normally an integral 
expression of the inverse square law) is to set up the Boundary Integral Equations (BIE) which 
lead to the moderately popular Boundary Element Method (BEM). In the forward version of 
the BEM, surfaces of a given geometry are replaced by a distribution of point singularities 
such as source / dipole of unknown strengths. The strengths of these singularities are obtained 
through the satisfaction of a given set of boundary conditions that can be Dirichlet, Neumann 
or of the Robin type. The numerical implementation requires considerable care pQ because it 
involves evaluation of singular (weak, strong and hyper) integrals. Some of the notable two- 
dimensional and three-dimensional approaches are pQ and El IU E] and the references in 
these papers. Despite a large body of literature, closed form analytic expressions for computing 
the effects of distributed singularities are rare [TJ E] and complicated to implement. Thus, for 
solving difficult but realistic problems involving, for example, sharp edges and corners or thin 
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elements, introduction of complicated mathematics and special formulations becomes a necessity 
[HJ[TU]. These drawbacks are some of the major reasons behind the relative unpopularity of the 
BEM despite its significant advantages over domain approaches such as the finite-difference and 
finite-element methods (FDM and FEM) while solving non-dissipative problems [TT]. It is well- 
understood that most of the difficulties in the available BEM solvers stem from the assumption 
of nodal concentration of singularities which leads to various mathematical difficulties and to the 
infamous numerical boundary layers [9], [12] . The Inverse Square Law Exact Solutions (ISLES) 
library, in contrast, is capable of truly modeling the effect of distributed singularities precisely 
and, thus, its application is not limited by the proximity of other singular surfaces or their 
curvature or their size and aspect ratio. The library consists of exact solutions for both potential 
and flux due to uniform distribution of singularity on flat rectangular and triangular elements. 
While the rectangular element can be of any arbitrary size [13\ I14j. the triangular element can 
be a right angled triangle of arbitrary size [TS]. Since any real geometry can be represented 
through elements of the above two types (or by the triangular type alone), this library can help 
in developing solvers capable of solving three-dimensional potential problems for any geometry. 
It may be noted here that any non-right-angled triangle can be easily decomposed in to two 
right-angled triangles. Thus, the right-angled triangles considered here, in fact, can take care 
of any three-dimensional geometry. Several difficulties were faced in developing the library 
which arose due to the various terms of the integrals and also from the approximate nature of 
computation in digital computers. In this paper, we have discussed these difficulties, solutions 
adopted at present and possible ways of future improvement. 

The classic benchmark problem of estimating the capacitance of a unit square plate raised 
to unit volt has been addressed using a solver based on ISLES, namely, the nearly exact BEM 
(neBEM) solver. Results obtained using neBEM have been compared with other precise results 
available in the literature. The comparison clearly indicates the excellent precision and efficiency 
achievable using ISLES and neBEM. In addition, we have also presented the variation of charge 
density close to the corner of the square plate. Usually, using BEM, it is difficult to obtain 
physically consistent results close to these geometric singularities. Wild variations in the magni- 
tude of the charge density has been observed with the change in the degree of discretization, the 
reason once again being associated with the nodal model of singularities [16]. In contrast, using 
neBEM, we have obtained very smooth variation close to the corner. Moreover, the magnitudes 
of the charge density have been found to be consistently converging to physically realistic values. 
These results clearly indicate that since the foundation expressions of the solver are exact, it 
is possible to find the potential and flux accurately in the complete physical domain, includ- 
ing the critical near-field domain using neBEM. In addition, since singularities are no longer 
assumed to be nodal and we have the exact expressions for potential and flux throughout the 
physical domain, the boundary conditions no longer need to be satisfied at special points such 
as the centroid of an element. Although consequences of this considerable advantage is still 
under study, it is expected that this feature will allow neBEM to yield accurate estimates for 
problems involving corners and edges that are very important in a large number of scientific and 
technological studies. 

It should be noted here that the exact expressions for triangular elements consist of a signifi- 
cantly larger number of mathematical operations than those for rectangular elements. Thus, for 
the solver, it is more economical if we use a mixed mesh of rectangular and triangular elements 
using rectangular elements as much as possible. However, in the present work, we have inten- 
tionally concentrated on the performance of the triangular elements and results shown here are 
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those obtained using only triangular elements. 



2 Exact Solutions 

The expressions for potential and flux at a point (X, Y, Z) in free space due to uniform source 
distributed on a rectangular flat surface having corners situated at {x\, z\) and (x%, Z2) has been 
presented, validated and used in [I3j [H] and, thus, is not being repeated here. 

Here, we present the exact expressions necessary to compute the potential and flux due to a 
right-angled triangular element of arbitrary size, as shown in FigHJ It may be noted here that 
the length in the X direction has been normalized, while that in the Z direction has been allowed 
to be of any arbitrary magnitude, zm- From the figure, it is easy to see that in order to find 
out the influence due to triangular element, we have imposed another restriction, namely, the 
necessity that the X and Z axes coincide with the perpendicular sides of the right-angled triangle. 
Both these restrictions are trivial and can be taken care of by carrying out suitable scaling and 
appropriate vector transformations. It may be noted here that closed- form expressions for the 
influence of rectangular and triangular elements having uniform singularity distributions have 
been previously presented in 018]. However, in these works, the expressions presented are quite 
complicated and difficult to implement. In [13] and in the present work, the expressions we have 
presented are lengthy, but completely straight-forward. As a result, the implementation issues 
of the present expressions, in terms of the development of the ISLES library and the neBEM 
solver are managed quite easily. 

m P 

Yi 




Figure 1: Right-angled triangular element with x- length 1 and an arbitrary z-length, zm', P is 
the point where the influence (potential and flux) is being computed 



It is easy to show that the influence (potential) at a point P(X, Y, Z) due to uniform source 
distributed on a right-angled triangular element as depicted in FigfT] can be represented as a 
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multiple of 

r 1 f z ( x ) dr dz 

in which we have assumed that x\ = 0, Z\ = 0, x 2 = 1 and z 2 = zm, as shown in the geometry 
of the triangular element. The closed-form expression for the potential has been obtained using 
symbolic integration [T7] which was subsequently simplified through substantial effort. It is found 
to be significantly more complicated in comparison to the expression for rectangular elements 
presented in [13] and can be written as 

$ = 

^ ( (z M Y 2 - XG){LP X + LM\ - LP 2 - LM 2 ) + i \Y\ (z M X + G){LP 1 - LM% - LP 2 + LM 2 ) 

+ iS 1 \Y\(tanh-H^±^) - tanh-\^-^) - tan h -\^±^) + 

2G yi + z M 2 Di2-E 1 ( D 21 -X + l \ 
+ h = g log ( /i Z 2 n — ) + 2Zlog — — )) + C (2) 

Vl + ZM VI + ZM P>2\ — Pi V - -A / 

where, 



D u = J(X- xi) 2 + + (Z - zi) 2 ; D 12 = J (X - Xl ) 2 + Y^ + (Z-z 2 



D 21 = ^(X- x 2 y + Yi + (Z- z x )A h = {X - x x ) \Y\ ; h = (X - x 2 ) \Y\ 
Si = sign(zi - Z); R 1 = Y 2 + (Z- z x ) 2 
E 1 = (X + zm 2 - z M Z); E 2 = (X-l- z M Z), 
G = z M (X - 1) + Z; H x = Y 2 + G(Z - z M ); H 2 = Y 2 + GZ 

1 , AH X + GD 12 )+i\Y\(E 1 -iz M D 12 ) 
LPl = G- iZM \Y\ l ° 9{ ^XTW\ } 



1 , AH 1 +GD 12 )-i\Y\(E l -iz M D l2 ) 



G + iz M \Y\ " v -X-i\Y\ 



1 . AH 2 +GD 21 )+i\Y\{E 2 -iz M D 21 ), 
LP2 = G- iZM \Y\ l ° 9{ l-X + i\Y\ } 



1 . AH 2 +GD 21 )-i\Y\(E 2 -iz M D 21 ] 
LM 2 = — ; ——log{- 



G + iz M \Y\ ^ l-X-i\Y\ ' 

and C denotes a constant of integration. 

Similarly, the flux components due to the above singularity distribution can also be repre- 
sented through closed-form expressions as shown below: 
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y dy 

^ ( {2z M Y){LP l + LMi - LP 2 - LM 2 ) + < |Y| (,Sn(y)G)(LP 1 - LMi - LP 2 + LM 2 ) 

(4) 

and, 



Fz = "^ = ( VTTlI? log( Vi + ^ 1 2-^ ) + log U 1 -x + i J ) + c 



(5) 



where Sn(Y) implies the sign of the Y-coordinate and C indicates constants of integrations. It 
is to be noted that the constants of different integrations are not the same. These expression 
are expected to be useful in the mathematical modeling of physical processes governed by the 
inverse square laws. Being exact and valid throughout the physical domain, they can be used to 
formulate versatile solvers to solve multi-scale multi-physics problems governed by the Laplace 
/ Poisson equations involving Dirichlet, Neumann or Robin boundary conditions. 



3 Development of the ISLES library 

Due to the tremendous popularity of the C language we have written the codes in the C pro- 
gramming language. However, it should be quite simple to translate the library to other popular 
languages such as FORTRAN or C++, since no special feature of the C language has been used 
to develop the codes. 



3.1 Validation of the exact expressions 

The expressions for the rectangular element have been validated in detail in [13] , Here, we present 
the results for triangular elements in fair detail. In FigEJ we have presented a comparison of 
potentials evaluated for a unit triangular element by using the exact expressions, as well as 
by using numerical quadrature of high accuracy. The two results are found to compare very 
well throughout. Please note that contours have been obtained on the plane of the element, 
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and thus, represents a rather critical situation. Similarly, Figj3] shows a comparison between 
the results obtained using closed-form expressions for flux and those obtained using numerical 
quadrature. The flux considered here is in the Y direction and is along a line beginning from 
(—2, —2, —2) and ending at (2, 2, 2). The comparison shows the commendable accuracy expected 

and |4(b)[ the surface plots of potential on the element 
plane (XZ plane) and Y-flux on the XY plane have been presented from which the expected 
significant increase in potential and sharp change in the flux value on the element is observed. 
Thus, by using a small fraction of computational resources in comparison to those consumed in 
numerical quadratures, ISLES can compute the exact value of potential and flux for singularities 
distributed on triangular elements. 



from closed form expressions. In Fig 4(a) 




Figure 2: Potential contours on a triangular element computed using exact expressions and by 
numerical quadrature 



3.2 Near-field performance 

In order to emphasize the accuracy of ISLES, we have considered the following severe situations 
in the near-field region in which it is observed that the quadratures can match the accuracy of 
ISLES only when a high degree of discretization is used. Please note that in these cases, the 
value of zm has been considered to be 10. In FigJ^we have presented the variation of potential 
along a line on the element surface running parallel to the Z-axis of the triangular element (see 
Fig{T]) and going through the centroid of the element. It is observed that results obtained using 
even a 100 x 100 quadrature is quite unacceptable. In fact, by zooming on to the image, it can 
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Figure 3: Comparison of flux (in the Y direction) as computed by ISLES and numerical quadra- 
ture along a diagonal line 



be found that only the maximum discretization yields results that match closely to the exact 
solution. It may be noted here that the potential is a relatively easier property to compute. 
The difficulty of achieving accurate flux estimates is illustrated in the two following figures. The 
variation of flux in the X-direction along the same line as used in FigJ5] has been presented in 
Figj6j Similarly, variation of Y-flux along a diagonal line (beginning at (-10,-10,-10) and ending 
at (10,10,10) and piercing the element at the centroid) has been presented in FigJTJ From these 
figures we see that the flux values obtained using the quadrature are always inaccurate even 
if the discretization is as high as 100 x 100. We also observe that the estimates are locally 
inaccurate despite the use of very high amount of discretization (200 x 200 or 500). Specifically, 
in the latter figure, even the highest discretization can not match the exact values at the peak, 
while in the former only the highest one can correctly emulate the sharp change in the flux 
value. It is also heartening to note that the values from the quadrature using higher amount of 
discretization consistently converge towards the ISLES values. 

3.3 Far field performance 

It is expected that beyond a certain distance, the effect of the singularity distribution can 
be considered to be the same as that of a centroidally concentrated singularity or a simple 
quadrature. The optimized amount of discretization to be used for the quadrature can be 
determined from a study of the speed of execution of each of the functions in the library and has 
been presented separately in a following sub-section. If we plan to replace the exact expressions 
by quadratures (in order to reduce the computational expenses, presumably) beyond a certain 
given distance, the quadrature should necessarily be efficient enough to justify the replacement. 
While standard but more elaborate algorithms similar to the fast multipole method (FMM) [TH] 
along with the GMRES [19J matrix solver can lead to further of computational efficiency, the 
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Figure 4: (a) Potential surface due to a triangular source distribution on the element plane, (b) 
Flux (in the Y direction) surface due to a triangular source distribution on the XY plane at Z=0 

simple approach as outlined above can help in reducing a fair amount of computational effort. 
In the following, we present the results of numerical experiments that help us in determining 
the far-field performance of the exact expressions and quadratures of various degrees that, in 
turn, help us in choosing the more efficient approach for a desired level of accuracy. 

In FiglH] we have presented potential values obtained using the exact approach, 100 x 100, 
10 x 10 and no discretization, i.e., the usual BEM approximation while using the zeroth order 
piecewise uniform charge density assumption. The potentials are computed along a diagonal line 
running from (-1000, -1000, -1000) to (1000, 1000, 1000) which pierces a triangular element of 
zm = 10. It can be seen that results obtained using the usual BEM approach yields inaccurate 
results as we move closer than distances of 10 units, while the 10 x 10 discretization yields 
acceptable results up to a distance of 1.0 unit. In order to visualize the errors incurred due 
to the use of quadratures, we have plotted Figj9] where the errors incurred (normalized with 
respect to the exact value) have been plotted. From this figure we can conclude that for the 
given diagonal line, the error due to the usual BEM approximation falls below 1% if the distance 
is larger than 20 units while for the simple 10 x 10 discretization, it is 2 units. It may be mentioned 
here that along the axes the error turns out to be significantly more [13] and the limits need 
to be effectively doubled to achieve the accuracy for all cases possible. Thus, for achieving 1% 
accuracy, the usual BEM is satisfactory only if the distance of the influenced point is five times 
the longer side of an element. Please note here that the error drops to 1 out of 10 6 as the 
distance becomes fifty times the longer side. Besides proving that the exact expressions work 
equally well in the near-field as well as the far-field, this fact justifies the usual BEM approach 
for much of the computational domain leading to substantial savings in computational expenses. 

The accuracy of the exact expressions used in the ISLES library is confirmed from the above 
comparisons. However, there are several other important issues related to the development of 
the library that are discussed below briefly. 

3.4 Evaluation of the component functions 

Many of the irrational and transcendental functions have domains and ranges in which they are 
defined. Moreover, they are often multiply defined in the complex domain; for example, there 
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Figure 5: Variation of potential along a centroidal line on the XZ plane parallel to the Z axis 
for a triangular element: comparison among values obtained using the exact expressions and 
numerical quadratures 



are an infinite number of complex values for the logarithm function. In such cases, only one 
principal value must be returned by the function. In general, such values cannot be chosen so 
as to make the range continuous and thus, lines in the domain called branch cuts need to be 
defined, which in turn define the discontinuities in the range. While evaluating expressions such 
as the ones displayed in eqns.([2]-E]) a number of such problems are expected to occur. However, 
when the expressions are analyzed at critical locations such as the corners and edges of the 
element, it is observed that the terms likely to create difficulties while evaluating potentials are 
either cancelled out or are themselves multiplied by zero. As a result, at these locations of likely 
geometric and mathematical singularities, the solution behaves nicely. However, the same is 
not true for the expressions related to the flux components. For these, we have to deal with 
branch-cut problems in relation to tanh' 1 and problems related to the evaluation of log(0). It 
should be noted here that these singularities associated to the edges and corners of the elements 
are of the weak type and it is expected that exact evaluation of these terms as well will be 
possible through further work. 

However, difficulties of a different nature crop up in these calculations which can be linked 
directly to the limitation of the computer itself, namely, round-off errors [TS] . These errors can 
lead to severe problems while handling multi-scale problems such as those described in [13]. A 
completely different approach is necessary to cope up with these difficulties, for example, the 
use of extended range arithmetic [20J, interval arithmetic [21J or the use of specialized libraries 
such as the CORE library of the Exact Geometric Computation (EGC) initiative [22]. In the 
present version of ISLES, a simple approach has been implemented which sets a lower limit to 
various distance values. Below this value, the distance is considered to be zero. Plan of future 
improvements in this regard has been kept at a high priority. 
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Figure 6: Variation of flux in the X direction along a line on the XZ plane parallel to the Z axis 
for a triangular element: comparison among values obtained using the exact expressions and 
numerical quadratures 



3.4.1 Algorithm 

As discussed above, there are possibilities of facing problems while using the exact expressions 
which may be due to the functions being evaluated or due to round-off errors leading to erro- 
neous results. Moreover, despite providing many checks during the computation there is finite 
possibility of ending up with a wrong value of a property indicated by its being Nan or inf 
or potential due to unit positive singularity strength turning out to be negative. In order to 
maintain the robustness of the library, we have tried to keep checks on the intermediate and final 
values during the course of the computation. When the results are found to be unsatisfactory, 
unphysical, we have re-estimated the results by using numerical quadrature and kept a track of 
the cause by raising a unique approximation flag which is specific for a problem. As a result, 
the steps for the calculation for a property can be written as follows: 

• Get the required inputs - geometry of the element and the position where the effect needs 
to be evaluated; Check whether the element size and distances are large enough so that 
the results do not suffer from round-off errors. 

• Check whether the location coincides with one of the special ones, such as corners or edges. 

• Evaluate the necessary expressions in accordance with the foregoing results. If necessary, 
consider each term in the expressions separately to sort out difficulties related to singular- 
ities, branch-cuts or round-off errors. Note that if the multiplier is zero, rest of the term 
does not need evaluation. 

• If direct evaluation of the expressions fail, raise a unique approximation flag specific to this 
problem and term and return the value of the property by using numerical quadrature. 
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Figure 7: Comparison of flux (in the Y direction) along a diagonal line piercing the triangular 
element at the centroid: comparison among values obtained using the exact expressions and 
numerical quadratures 



• Compute all the terms and find the final value, Check whether the final value is a number 
and physically meaningful. If not satisfied, recompute the result using numerical quadra- 
ture and raise the relevant approximation flag. 

3.5 Speed of execution 

The time taken to compute the potential and flux is an important parameter related to the 
overall computational efficiency of the codes. This is true despite the fact that, in a typical 
simulation, the time taken to solve the system of algebraic equation is far greater than the time 
taken to build the influence coefficient matrix and post-processing. Moreover, the amount of 
time taken to solve the system of equations tend to increase at a greater rate than the time 
taken to complete the other two. It should be mentioned here that the time taken in each of 
these steps can vary to a significant amount depending on the algorithm of the solver. In the 
present case, the system of equations has been solved using lower upper decomposition using 
the well known Crout's partial pivoting. Although this method is known to be very rugged 
and accurate, it is not efficient as far as number of arithmetic operations, and thus, time is 
concerned. It is also possible to reduce the time taken to pre-process (generation of mesh and 
creation of influence matrices), solve the system of algebraic equations and that for post-process 
(computation of potential and flux at the required locations) can be significantly reduced by 
adopting faster algorithms, including those involving parallelization. 

In order to optimize the time taken to generate the influence coefficient matrix and that to 
carry out the post-processing, we carried out a small numerical study to determine the amount 
of time taken to complete the various functions being used in ISLES, especially those being used 
to evaluate the exact expressions and those being used to carry out the quadratures. The results 
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Figure 8: Potential along a diagonal through the triangular element computed using exact, 
100 x 100, 10 x 10 and usual BEM approach 



of the study (which was carried out using the linux system command gprof) has been presented 
in the following Tabl^H 



Table 1: Time taken to evaluate exact expression and various quadratures 



Method 


Exact 


Usual BEM 


10 x 10 


100 x 100 


500 x 500 


Time 


0.8 /is 


25 ns 


1 /is 


200 lis 


5 ms 



Please note that the numbers presented in this table are representative and are likely to 
have statistical fluctuations. However, despite the fluctuations, it may be safely concluded that 
a quadrature having only 10 x 10 discretization is already consuming time that is comparable 
to that needed exact evaluation. Thus, the exact expressions, despite their complexity, are 
extremely efficient in the near-field which can be considered at least as large as 0.5 times the 
larger side of a triangular element (please refer to Figl9]). In making this statement, we have as- 
sumed that the required accuracy for generating the influence coefficient matrix and subsequent 
potential and flux calculations is 1%. This may not be acceptable at all under many practical 
circumstances, in which case the near-field would imply a larger volume. 

3.6 Salient features of ISLES 

Development of usual BEM solvers are dependent on the two following assumptions: 

• While computing the influences of the singularities, the singularities are modeled by a sum 
of known basis functions with constant unknown coefficients. For example, in the constant 
element approach, the singularities are assumed to be concentrated at the centroid of the 
element, except for special cases such as self influence. This becomes necessary because 
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Figure 9: Error along a diagonal through the triangular element computed using 100 x 100, 
10 x 10 and usual BEM approach 



closed form expressions for the influences are not, in general, available for surface elements. 
An approximate and computationally rather expensive way of circumventing this limitation 
is to use numerical integration over each element or to use linear or higher order basis 
functions. 

• The strengths of the singularities are solved depending upon the boundary conditions, 
which, in turn, are modeled by the shape functions. For example, in the constant element 
approach, it is assumed that it is sufficient to satisfy the boundary conditions at the 
centroids of the elements. In this approach, the position of the singularity and the point 
where the boundary condition is satisfied for a given element usually matches and is called 
the collocation point. 

The first (and possibly, the more damaging) approximation for BEM solvers can be relaxed by 
using ISLES and can be restated as, 

• The singularities distributed on the boundary elements are assumed to be uniform on a 
particular element. The strength of the singularity may change from element to element. 

This improvement turns out to be very significant as demonstrated in the following section and 
some of our other studies involving microelectromechanical systems (MEMS) and gas detectors 
for nuclear applications [EHH]. Some of the advantages of using ISLES are itemized below: 

• For a given level of discretization, the estimates are more accurate, 

• Effective efficiency of the solver improves, as a result, 

• Large variation of length-scales, aspect ratios can be tackled, 
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• Thinness of members or nearness of surfaces does not pose any problem, 

• Curvature has no detrimental effect on the solution, 

• The boundary condition can be satisfied anywhere on the elements, i.e., points other than 
the centroidal points can be easily used, if necessary (for a corner problem, may be), 

• The same formulation, library and solver is expected to work in majority of physical 
situations. As a result, the necessity for specialized formulations of BEM can be greatly 
minimized. 

4 Capacitance of a unit square plate - a classic benchmark prob- 
lem 

Using the neBEM solver, we have computed the capacitance of a unit square conducting plate 
raised to a unit volt. This problem is still considered to be one of the major unsolved problems of 
electrostatic theory [231 IS, 26, 16J and no analytical solution for this problem has been obtained 
so far. The capacitance value estimated by the present method has been compared with very 
accurate results available in the literature (using BEM and other methods). The results obtained 
using the neBEM solver is found to be among the most accurate ones available till date as shown 
in TableEl Please note that we have not invoked symmetry or used extrapolation techniques to 
arrive at our result presented in the table. 



Table 2: Comparison of capacitance values 



Reference 


Method 


Capacitance (pF) / 4 ireo 


[23j 


Surface Charge 


0.3607 


[24j 


Surface Charge 


0.362 


[25j 


Surface Charge 


0.367 


i 


Refined Surface Charge 


0.3667892 ± 1.1 x 10~ 6 




and Extrapolation 




[26] 


Refined Boundary Element 
and Extrapolation 


0.3667874 ±lx 10~ 7 


m 


Numerical Path Integration 


0.36684 


[161 


Random Walk 


0.36 ± 0.01 


This work 


neBEM 


0.3660587 



Finally, we consider the corner problem related to the electrostatics of the above conducting 
plate. Problems of this nature are considered to be challenging for any numerical tool and 
especially so for the BEM approach. The inadequacy of the BEM approach, especially in 
solving the present problem, has been mentioned even quite recently [16j in which it has been 
correctly mentioned that since the method can not extend its mathematical model to include 
the edges and corners in reality, it is unlikely that it will ever succeed in modeling the edge 
/ corner singularities correctly. As a result, with change in discretization, the properties near 
these geometric singularities are expected to oscillate significantly leading to erroneous results. 
However, as discussed above, the neBEM does extend its singularities distributed on the surface 
elements right till an edge or a corner. Moreover, using neBEM, it is also possible to satisfy the 
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boundary conditions (both potential and flux) as close to the edge / corner as is required. In 
fact, it should be possible to specify the potentials right at the edge / corner. 

In the following study, we have presented estimates of charge density very close the flat plate 
corner as obtained using neBEM. Please note that the boundary conditions have been satisfied 
at the centroids of each element although we plan to carry out detailed studies of changing the 
position of these points, especially in relation to problems involving edges / corners. In Fig llOl 
charge densities very close to the corner of the flat plate estimated by neBEM using various 
amounts of discretization have been presented. It can be seen that each curve follows the same 
general trend, does not suffer from any oscillation and seems to be converging to a single curve. 
This is true despite the fact that there has been almost an order of magnitude variation in the 
element lengths. 

Finally, in Fig llll we present a least-square fitted straight line matching the charge density 
as obtained the highest discretization in this study. It is found that the slope of the straight 
line is 0.713567, which compares very well with both old and recent estimates of 0.7034 |28} 129], 
This is despite the fact that here we have used a relatively coarse discretization near the corner. 
It should be mentioned that none of the earlier references cited here used the BEM approach. 
While the former used a singular perturbation technique, the latter used a diffusion based Monte- 
Carlo method. Thus, it is extremely encouraging to note that using the neBEM approach, we 
have been able to match the accuracy of these sophisticated techniques. 




Figure 10: Corner charge density estimated by neBEM using various sizes of triangular elements 



5 Conclusion 

An efficient and robust library for solving potential problems in a large variety of science and 
engineering problems has been developed. Exact closed-form expressions used to develop ISLES 
have been validated throughout the physical domain (including the critical near-field region) 
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Figure 11: Variation of charge density with increasing distance from the corner of the unit square 
plate and a least-square fitted straight line: slope of the fitted line is 0.713567 



by comparing these results with results obtained using numerical quadrature of high accuracy. 
Algorithmic aspects of this development have also been touched upon. A classic benchmark 
problem of electrostatics has been successfully simulated to very high precision. Charge density 
values at critical geometric locations like corners have been found to be numerically stable and 
physically acceptable. Several advantages over usual BEM solvers and other specialized BEM 
solvers have been briefly mentioned. Work is under way to make the code more robust and 
efficient through the implementation of more efficient algorithms and parallelization. 
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